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ABSTRACT 

Over  recent  years,  a  variety  of  shock-capturing  schemes  have  been  developed  for  the 
Euler  equations  of  gas  dynamics.  During  this  period,  it  has  emerged  that  one  of  the  more 
successful  strategies  is  to  follow  Godunov’s  lead  and  utilize  a  nonlinear  building  block  known 
as  a  Riemann  problem.  Now,  although  Rietnaun  solver  technology  is  often  thought  of  as 
being  mature,  there  are  in  fact  several  circumstances  for  which  Godunov-type  schemes  are 
found  wanting.  Indeed,  one  inherent  deficiency  is  so  severe  that  if  left  unaddressed,  it  could 
preclude  such  schemes  from  being  used  to  capture  detonation  fronts  in  simulations  of  complex 
flow  phenomena.  In  this  paper,  we  highlight  this  particular  deficiency  along  with  some  other 
little  known  weaknesses  of  Godunov-type  schemes,  and  we  outline  one  strategy  that  we  have 
used  to  good  effect  in  order  to  produce  reliable  high  resolution  simulations  of  both  reactive 
and  nonreactive  shock  wave  phenomena.  In  particular,  we  present  results  for  simulations  of 
so-called  galloping  instabilities  and  detonation  cell  phenomena. 
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1.  Introduction 


Following  the  rise  of  computational  fluid  dynamics,  numerical  simulations  of  shock  wave 
phenomena  are  now  commonplace.  Such  simulations  are  attractive  as  replacements  for  ex¬ 
periments  which  are  either  difficult,  dangerous,  or  expensive,  and  can  be  done  for  problems 
which  are  not  amenable  to  analytical  methods.  However,  because  of  the  disparate  scales 
involved,  simulations  are  all  too  often  under- resolved  and  so  are  of  limited  use.  Indeed, 
despite  the  plethora  of  numerical  schemes  that  have  been  developed,  only  the  Godunov-type 
methods  have  been  shown  to  produce,  genuinely,  high  fidelity  simulations  of  complex  shock 
wave  phenomena  (Berger  and  Colella,  1989;  Quirk,  1992  a).  Consequently,  following  a  survey 
of  the  literature,  Godunov-type  methods  are  likely  to  be  picked  up  by  people  who  are  not 
algorithm  developers  as  tools  with  which  to  simulate  real  problems.  Unfortunately,  despite 
their  undoubted  strengths,  Godunov-type  schemes  have  certain  inherent  weaknesses  that 
can,  occasionally,  result  in  simulations  failing  catastrophically.  Amongst  the  cognoscenti, 
such  weaknesses,  whilst  not  always  fully  understood,  can  be  overcome.  However,  the  non¬ 
expert  is  severely  disadvantaged  by  the  fact  that  the  various  failings  and  their  associated 
fixes  go  largely  unreported.  Recently  we  have  attempted  to  redress  this  unfortunate  state  of 
affairs  (Quirk,  1992  b ),  and  here  we  present  an  abridged  version  of  this  work  which  is  targeted 
directly  at  people  within  the  combustion  community  who  would  like  to  take  advantage  of 
modern  shock-capturing  schemes  for  simulating  detonation  flows. 

The  rest  of  this  paper  is  as  follows.  In  Section  2  we  present  a  brief  outline  of  Godunov- 
type  schemes.  This  is  followed  by  a  section  containing  some  specific  examples  of  how  different 
schemes  can  fail.  In  Section  4  we  describe  the  strategy  that  we  use  to  improve  the  robustness 
of  our  flow  solvers,  following  which,  we  present  results  for  simulations  of  galloping  instabilities 
and  detonation  cell  phenomena.  Finally,  in  Section  6  we  present  some  conclusions  that  we 
have  drawn  from  this  work. 

2.  Outline  of  Godunov-type  Schemes 

Many  expositions  of  Godunov’s  method  and  its  descendants  appear  in  the  literature 
(Holt,  1984;  Roe,  1986);  here  we  simply  want  to  present  the  general  gist  of  such  schemes  in 
order  to  orientate  the  reader  for  the  material  which  follows  in  the  remaining  sections  of  this 
paper. 

With  reference  to  Figure  1,  a  Godunov-type  scheme  may  be  viewed  as  follows.  The 
scheme  works  with  a  low-order  projection  of  the  flow  solution;  each  mesh  cell  contains  a 
cell-averaged  value  for  the  true  solution  over  the  cell.  Thus,  the  numerical  representation 
closely  approximates  the  true  solution  near  discontinuities,  and  regions  of  smooth  flow  are 
reasonably  well  approximated  by  a  series  of  step  functions.  This  discrete  system  is  inte¬ 
grated  by  first  reconstructing  the  flow  solution  within  each  cell.  This  step  is  effectively  an 
extrapolation  process  for  finding  the  flow  states  at  the  edges  of  mesh  cells  given  values  at 
the  centres  of  cells.  Note  that,  in  general,  the  reconstructed  solution  whilst  smoother  than 
the  projected  solution  will  still  be  discontinuous  at  cell  interfaces.  Next  there  comes  an 
evolution  step.  A  Riemann  problem  is  solved  for  each  cell  interface  using  the  reconstructed 
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Cell-averaged  projection  of  the 
flow  solution. 


Pierewise-linear  reconstruction 
of  the  flow  solution. 


Typical  Riemann  solution 
showing  three  centred  waves. 


Figure  1:  A  Godunov-type  scheme  consists  of  three  basic  steps:  reconstruction— + evolution 
— > projection . 


states  on  either  side  of  the  interface  as  input  data.  Recall  that  the  Riemann  problem  for 
any  set  of  conservation  laws  arises,  if  initial  data  are  prescribed  as  two  semi-infinite  states 
(W  =  Wi  for  x  <  Q,  W  =  W r  for  x  >  0).  The  solution  to  a  Riemann  problem  consists  of  a 
number  of  centred  waves.  For  the  one-dimensional  Euler  equations  of  gas  dynamics  there 
are  three  waves.  The  inner  wave  is  a  contact  discontinuity  separating  states  at  different 
temperatures,  and  each  of  the  two  outer  waves  may  be  either  a  shock  wave  or  an  expansion 
wave.  Finally,  the  different  solutions  to  the  separate  Riemann  problems  are  averaged  so  as 
to  find  a  cell-centred  projection  of  the  flow  solution  at  some  new  time  level.  If  repeated, 
the  sequence  reconstruction— ►  evolution  — ► projection  results  in  an  accurate  and  well-behaved 
scheme  for  simulating  the  propagation  of  shock  waves. 

Whilst  the  overall  strategy  of  a  Godunov-type  scheme  is  largely  clear-cut2,  the  same  can¬ 
not  be  said  of  its  individual  components.  For  example,  more  often  than  not,  the  Riemann 
problem  for  a  system  of  conservation  laws  does  not  have  an  analytic  solution  and  is  there¬ 
fore  expensive  to  compute.  Consequently,  many  workers  prefer  to  compute  an  approximate 
uohd'on  to  the  Riemann  problem  which  embodies  the  spirit  of  the  exact  solution  but  which 
is  cheaper  to  compute  (Roc,  19^8;  Einfeldt,  1988;  Toro,  1991).  Indeed,  the  design  of  approx¬ 
imate  Riemann  solvers  has  become  something  of  an  industry  in  its  own  right,  and  there  is 
considerable  debate  concerning  the  relative  merits  of  different  solvers.  Moreover,  there  are 

2 Here  we  have  describe'}  the  so  r piled  Ml'SCI.  approach  f. ;  pro  facing  a  high-order  Godunov  method. 
An  alternative  methodology  is  followed  by  flux-limited  schemes  where  the  reconstruction  step  is  replaced  by 
a  procedure  which  post- processes  the  Riemann  solution. 
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many  circumstances  for  which  certain  approximate  solvers  actually  prove  to  be  more  reliable 
than  the  exact  solver.  Therefore,  in  general,  it  is  far  from  obvious  which  particular  Riemann 
solver  is  best  suited  to  a  given  application. 

The  reconstruction  step  is  similarly  open  to  different  interpretation.  Firstly,  there  is  a 
free  choice  as  to  which  variables  are  reconstructed,  and  which  quantities  are  derived  from 
the  reconstructed  variables.  For  example,  although  shock  capturing  schemes  invariably  work 
w’ith  the  conserved  variables,  experience  shows,  at  least  for  the  Euler  equations,  that  better 
results  are  obtained  if  the  reconstruction  is  performed  using  primitive  variables.  Secondly, 
there  is  also  a  choice  as  to  the  order  of  accuracy  of  the  reconstruction.  As  is  common  practice, 
we  employ  piecewise-linear  slopes  (Quirk, 1992  a);  however,  C'olella  and  Woodward  (1984) 
use  piecewise-parabolas  to  good  effect  in  their  PPM  scheme,  and  more  recently  Harten  rt 
al.  (1987)  have  introduced  the  idea  of  ENO  schemes  where  the  reconstruction  step  may  be 
carried  out  to  some  arbitrary  order  of  accuracy. 

It  is  worth  noting,  however,  that  to  a  large  extent  the  quality  of  results  produced  will 
depend  on  how  well  the  flow  solution  is  reconstructed  near  discontinuities  and  so  the  notional 
order  of  accuracy  for  some  reconstruction  process  is  not  necessarily  a  reliable  indicator  of  its 
actual  performance.  Typically,  limiter  functions  (Roe,  1985)  are  employed  so  as  to  ensure 
that  the  reconstruction  process  does  not  introduce  unwanted  overshoots  near  extrema,  and 
it  is  the  properties  of  the  chosen  limiter  function  that  largely  dictate  the  quality  of  the  flow 
solution.  But  a  limiter  function  is  only  as  good  as  the  data  upon  which  it  acts.  For  example, 
although  we  employ  a  piecewise-linear  reconstruction  which  results  in  a  scheme  that  is  nom¬ 
inally  only  second-order  accurate,  the  slopes  are  derived  from  the  Riemann  solutions  which 
are  computed  using  the  cell-centred  states  as  input  data.  For  very  high  resolution  simula¬ 
tions  which  involve  thousands  of  time  steps,  this  strategy  gives  markedly  better  results  for 
weak  discontinuities  such  as  contact  surfaces  than  does  the  third-order  reconstruction  pro¬ 
cess  proposed  by  Anderson  et  al.(  1985)  which  does  not  utilize  the  same  level  of  characteristic 
information.  While  the  differences  are  much  less  marked  for  problems  that  contain  just  a 
few  time  steps,  they  are  still  appreciable  (see  Figure  2). 

Everything  considered,  given  the  basic  framework  of  a  Godunov-type  scheme,  it  is  possible 
to  construct  many  different  variations  on  a  theme.  Unfortunately,  there  is  no  “correct"  way  of 
doing  things,  for  the  theory  which  underpins  this  class  of  scheme  does  not  always  discriminate 
between  the  different  options  that  are  available.  In  practice,  seemingly  innocuous  details  of 
the  reconstruction  process  can  have  a  large  bearing  on  the  quality  of  the  results  produced  for 
the  sorts  of  very  detailed  simulations  that  will  be  necessary  in  order  to  unravel  the  dynamics 
of  detonation  phenomena.  Thankfully,  the  large  scale  flow  features  are  normally  insensitive 
to  such  changes,  and  so  one’s  faith  in  Godunov-type  schemes  is  not  unduly  undermined  by 
the  uncertainties  in  the  precise  details  of  the  method. 

Before  proceeding  to  the  next  section  which  exposes  some  of  the  failings  of  Riemann 
solve;.,,  so  as  not  to  appear  to  paint  an  overly  pessimistic  view  of  the  capabilities  of  Godunov- 
type  schemes,  it  is  worthwhile  showing  what  can  be  achieved  if  due  care  is  taken.  Figure  ;j 
shows  a  snapshot  taken  from  the  simulation  of  a  planar  shock  wave  diffracting  around  a  90° 
corner  (Quirk,  1992  b).  This  picture  is  similar  to  a  Schlieren  image  in  that  the  different 
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Figure  2:  The  internal  energy  profiles  computed  for  Sod's  problem  using:  (a)  2,l,i  order 
characteristic  MUSCL  (Quirk,  1992).,  (b)  'Vd  order  non-characteristic  MUSCL  (Anderson  ft 
al,  1985). 

shades  of  grey  depict  the  magnitude  of  the  gradient  of  the  density  field  (the  darker  the 
shade,  the  larger  the  gradient),  and  it  clearly  shows  all  the  salient  flow  features  as  identified 
experimentally  by  Bazeuhova  et.  al.  (1984). 

Briefly,  with  reference  to  Figure  4,  the  diffraction  of  the  incident  shock  wave  (AC)  around 
the  corner  gives  rise  to  an  expansion  fan  which  emanates  from  (O).  The  shape  of  this  fan's 
lead  characteristic  (AQO)  indicates  that  the  flow  upstream  of  the  incident  shock  is  super¬ 
sonic.  The  expansion  fan  interacts  with  the  incident  shock  to  form  the  disturbed  shock  front 
(AL)MN).  The  incident  shock  is  sufficiently  strong  that  this  disturbed  front  is  kinked;  a  Mach 
reflection  at  the  wall  gives  rise  to  a  triple  point  at  (M).  A  contact  surface  (ALO)  marks  the 
boundary  between  fluid  that  has  been  induced  into  motion  by  the  incident  shock  and  fluid 
which  has  been  processed  by  the  disturbed  shock  front.  Note  that  the  flow  is  separated  from 
the  wall  at  a  point  slightly  downstream  of  the  apex  of  the  corner,  and  a  slipstream  (OS) 
separates  the  expanded  flow  from  this  region  of  almost  stationary  gas.  The  free  end  of  this 
slipstream  rolls  up  into  a  vortex.  A  secondary  shock  wave  (KTS)  matches  the  pressure  of  the 
flow  accelerated  by  the  expansion  fan  to  that  of  the  decelerated  flow  behind  the  disturbed 
part  of  the  incident  wave.  This  secondary  shock  is  kinked  as  a  result  of  its  interaction  with 
the  slip  stream  (OS)..  A  secondary  contact  surface  (TL)  begins  at  the  point  of  intersection 
of  the  secondary  shock  wave  with  the  weak  shock  (OT)  which  terminates  the  expansion 
fan.  Lastly,  a  shock  wave  (PB)  is  present  so  as  to  decelerate  the  reversed  flow  within  the 
separated  region  as  it  approaches  the  point  of  diffraction. 
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3.  Some  Failings  of  Riemann  Solvers 

YVe  now  present  some  examples  where  certain  Riemann  solvers  are  known  to  give  unreli¬ 
able  results  for  the  Euler  equations  of  gas  dynamics.  Whilst  our  catalogue  of  failings  is  not 
exhaustive,  it  should  alert  the  reader  to  the  types  of  problems  that  they  might  encounter 
when  using  Clodunov-type  methods.  It  is  important  to  realise  that  no  one  failing  afflicts 
all  Riemann  solvers.  Conversely,  however,  it  would  appear  that  no  one  Riemann  solver  is 
completely  free  of  defects. 

With  reference  to  Figure  5,  Example  (a)  shows  how  spurious  post-shock  oscillations  can 
occur  whenever  a  shock  wave  moves  very  slowly  during  the  course  of  a  simulation.  Here, 
a  one-dimensional  shock  wave  is  moving  from  left  to  right  (for  a  Courant  number  of  one. 
it  takes  approximately  50  time  steps  for  the  shock  to  traverse  a  single  cell).  As  the  shock 
moves  relative  to  the  mesh,  so  the  smeared  numerical  shock  profile  inevitably  changes  shape. 
Unfortunately,  for  many  Riemann  solvers  including  the  exact  solver,  the  points  within  the 
smeared  profile  do  not  lie  on  the  Hugoniot  curve  connecting  the  pre-shock  state  to  the  post¬ 
shock  state.  Thus  these  two  states  cannot  be  connected  by  a  single  family  of  acoustic  waves. 
For  such  schemes,  whenever  the  shock  profile  changes,  so  the  supposedly  passive  wave  fields 
are  activated  (in  this  case,  the  u  and  u  —  a  wave  fields),  thus  giving  rise  to  low  frequency 
oscillations.  A  thorough  description  of  this  particular  failing  has  been  given  by  Roberts 
(1990).  Note  that  these  oscillations  are  not  the  same  as  the  high-frequency,  post-shock 
oscillations  that  afflict  finite-difference  shock-capturing  schemes. 

Example  (b)  in  Figure  5  is  taken  from  a  simulation  of  double  Mach  reflection  which 
was  performed  using  Roe’s  scheme;  here  we  show  a  snapshot  of  the  pressure  contours  at 
one  instant  in  time.  The  Mach  stem  is  inexplicably  kinked  giving  rise  to  a  spurious  triple 
point  (S).  It  should  be  noted  that  this  kinking  is  not  related  to  the  slight  bulging  that  is 
often  observed  experimentally  for  this  type  of  shock  reflection  problem.  Such  bulging  arises 
because  the  contact  discontinuity  emanating  from  the  primary  triple  point  is  deflected  by  the 
reflecting  surface  resulting  in  a  strong  wall  jet  which  effectively  pushes  out  the  base  of  the 
Mach  stem.  The  mechanism  behind  this  failing  is  not  fully  understood,  but  it  would  appear 
to  arise  from  the  fact  that  the  Mach  stern  is  closely  aligned  with  the  computational  grid. 
Consequently,  little  or  no  dissipation  is  provided  by  Roe’s  scheme,  in  a  direction  parallel  to 
the  stem,  to  control  the  kinking. 

Generally  speaking,  the  dissipation  mechanism  provided  by  many  Riemann  solvers  proves 
to  be  inadequate  whenever  a  strong  shock  wave  is  aligned  with  the  computational  grid. 
Example  (c)  shows  the  so-called  carbuncle  phenomena  (Peery  and  Imlay,  1988)  where  some 
schemes  fail  to  produce  a  realistic  bow  shock  for  a  blunt  body  placed  in  a  high  Mach  number 
flow  (here  we  have  plotted  density  contours).  Note  that  along  the  stagnation  line  the  bow 
shock  is  more  or  less  aligned  with  the  body-fitted  grid  used  for  the  calculation,  and  so  very 
little  dissipation  is  added  normal  to  the  stagnation  line  in  the  vicinity  of  the  bow  shock;  a 
small  amount  of  auxiliary  dissipation  applied  in  this  region  is  generally  suffic’ent  to  suppress 
the  carbuncle. 

Example  (d)  in  Figure  5  shows  a  particularly  insidious  failing  that  can  occur  when  a 
strong  shock  is  aligned  with  the  grid.  Here  we  have  plotted  a  single  snapshot  of  the  density 


(d) 


Figure  5:  Some  Riemann  solver  failings:  (a)  Slowly  moving  shocks.,  (b)  Kinked  Mach  stems., 
(c)  The  Carbuncle  phenomena.,  (d)  Odd- Even  decoupling. 
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contours  taken  from  the  simulation  of  an  initially  planar  shock  wave  which  is  propagating 
down  a  duct,  from  left  to  right.  A  nominally  uniform  Cartesian  grid  was  used  for  the 
calculation;  the  grid  centre-line  carried  a  small  saw-tooth  perturbation.  This  perturbation 
acts  like  a  forcing  function  which  causes  odd  even  decoupling  to  occur,  along  the  length  of 
the  shock,  in  both  the  density  and  pressure  fields.  Interestingly,  within  the  body  of  the 
shock,  the  decoupling  of  the  density  field  is  out  of  phase  with  that  of  the  pressure  field.  An 
attempt  to  explain  this  failing  has  been  given  by  Quirk  (1992  />),  but  certain  details  of  the 
mechanism  remain  unclear.  It  is  clear,  however,  that  this  numerical  instability  only  occurs 
for  strong  shocks,  and  that  it  takes  a  tong  time  to  develop  (here,  the  shock  had  propagated 
approximately  thirty  channel  widths  before  the  instability  first  became  apparent ).  Moreover, 
as  will  be  shown  in  Section  5,  this  particular  instability  lias  quite  grave  consequences  for  the 
simulation  of  two-dimensional  detonation  fronts;  waves  produced  by  the  numerical  instability 
can  interfere  with  the  genuine  transverse  waves  associated  with  the  propagation  of  the  front. 

For  many  problems,  the  nature  of  the  flow  solution  is  known  a  priori,  and  so  it  is  fairly 
obvious  whenever  a  simulation  gives  anomalous  results.  But  often  no  such  safety  net  exists,  in 
which  case  it  becomes  very  difficult  to  determine  the  fidelity  of  a  numerical  simulation.  One 
tactic  that  we  routinely  employ  is  to  run  a  simulation  two  or  more  times,  each  time  varying 
the  elements  of  the  flow  solver.  For  example,  we  may  change  our  choice  of  Riemann  solvers, 
or  we  may  change  our  choice  of  variables  for  the  reconstruction  process.  Admittedly,  the  fact 
that  two  such  simulations  give  similar  behaviour  is  no  guarantee  that  the  results  are  correct, 
but  it  is  useful  for  determining  which  features  of  the  solution  are  likely  to  be  numerical 
artifacts.  Lastly,  it  should  he  noted  that  many  of  the  failings  associated  with  Godunov-type 
methods  could  be  circumvented  if  strong  shocks  were  fitted  rather  than  captured.  However, 
given  the  complexity  of  a  general  purpose  shock-fitting  scheme,  this  option  is  not  likely  to 
appeal  to  the  worker  who  is  merely  using  a  Godunov- type  scheme  as  a  tool. 


4.  An  Adaptive  Riemann  Solver 

Having  exposed  some  of  the  weaknesses  of  Riemann  solvers,  we  now  present  a  simple 
strategy,  that  we  have  found  useful,  for  improving  the  all-round  robustness  of  Godunov-type 
schemes.  In  essence,  we  select  the  precise  flavour  of  upwinding  to  match  the  local  flow  data 
such  that  a  particular  Riemann  solver  is  only  employed  in  those  situations  where  it  is  known 
to  give  reliable  results.  By  recognizing  the  limitations  of  any  one  solver  it  is  possible  to  reap 
its  advantages  without  suffering  its  attendant  failings. 

Our  synergetic  strategy  has  a  number  of  attractions,  not  least  of  which  is  that  some 
favoured  solver  need  not  be  jettisoned  simply  because  it,  occasionally,  fails.  However,  it 
does  introduce  the  difficulty  of  how  to  decide  when  to  use  one  Riemann  solver  in  preference 
to  another.  But  it  has  been  our  experience  that  this  added  difficulty  is  not  particularly 
bothersome,  for  we  tend  to  combine  a  single  high  resolution  Riemann  solver  with  just  one 
or  two  other  solvers  that  prove  more  reliable  under  conditions  which  are  fairly  well-defined, 
and  so  a  set  of  ad  hoc  switching  functions  suffice.  For  example,  some  of  the  worst  failings 


of  Riemann  solvers  occur  in  the  vicinity  of  strong  shock  waves.  To  overcome  such  failings 
we  use  the  HLLE  scheme  (Einfeldt,  1988).  Now  it  makes  little  sense  to  chop  and  change 
the  choice  of  Riemann  solver  used  along  the  length  of  a  shock  wave,  since  to  do  so  would 
inevitably  perturb  a  planar  shock  from.  Hence,  we  apply  this  particular  Riemann  solver 
throughout  the  immediate  vicinity  of  a  strong  shock.  1  bus  the  HLLE  switching  function 
need  only  locate  the  position  of  a  shock  wave,  but  such  functions  already  exist  in  the  guise 
of  mesh  refinement,  monitor  functions. 

A  simple  test  that  identifies  those  cell  interfaces  which  are  in  the  vicinity  of  a  strong 
shock  is  to  check  whether  or  not 


\Pr  ~  Pi  I 


’nin(pt,pT) 


where  a  is  some  threshold  parameter  which  is  problem  dependent  and  pr  and  pi  refer  to 
the  pressures  which  act  on  the  interface.  If  this  condition  is  met,  the  two  cells  separated  by 
the  interface  are  flagged  as  lying  within  a  strong  shock.  Then,  when  it  comes  to  computing 
cell-interface  fluxes,  if  the  cells  either  side  of  an  interface  have  both  been  flagged  as  lying 
within  a  strong  shock,  the  flux  is  computed  using  the  HLLE  solver.  Note  that  since  numerical 
shocks  are  invariably  smeared  over  several  mesh  cells,  it  is  worth  locating  shocks  using  a 
projection  of  the  flow  solution  on  a  grid  which  is  coarser  than  that  used  for  the  calculation. 
On  such  a  grid  a  shock  will  appear  much  less  smeared,  and  so  the  left-hand  side  of  the  above 
switching  function  will  be  a  fair  indication  of  its  strength.  Once  a  set  of  cells  have  been 
flagged  on  this  coarse  mesh,  the  flags  may  be  prolongated  to  the  actual  computational  mesh 
so  as  to  find  those  cells  which  lie  in  the  vicinity  of  a  shock  wave. 

Figure  6  shows  how  the  HLLE  solver  may  be  used  to  correct  the  tendency  of  Roe's  scheme 
to  produce  kinked  Mach  stems,  c.f.  Figure  5  (b).  For  this  calculation  the  HLLE  switching 
function  was  tuned  such  that  it  would  only  be  activated  by  the  incident  shock,  and  the 
principal  Mach  stem;  the  threshold  parameter  n  was  simply  set  to  half  the  strength  of  the 
incident  shock  wave  as  given  by  the  left-hand  side  of  Equation  1.  Note  that  apart  from  the 
region  near  the  Mach  stem,  these  new  results  are  very  similar  to  the  old  ones.  This  shows 
that  the  HLLE  scheme  has  had  no  adverse  affect  on  the  resolution  of  Roe's  scheme. 

Having  presented  the  gist  of  our  strategy,  we  see  little  point  in  trying  to  sell  a  particular 
combination  of  solvers.  Starting  with  some  high  resolution  Riemann  solver,  whose  c.ioice  will 
inevitably  be  a  matter  of  personal  taste,  the  correct  combination  of  solvers  will  depend  both 
on  that  schemes  weaknesses  and  on  the  specific  application  in  hand.  In  turn,  the  combination 
of  Riemann  solvers  will  dictate  the  choice  of  switching  functions. 
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(a) 


(b) 


Figure  6:  The  HLLE  scheme  can  be  uied  to  circumvent  the  tendency  of  Roe's  method  to 
produce  kinked  Mach  stems:  (a)  Pressure  Contours,  (b)  HLLE  switching  function. 


5.  Galloping  Instabilites  and  Detonation  Cells 

Almost  all  of  our  experience  with  Godunov-type  schemes  has  been  gained  from  calcula¬ 
tions  of  nonreactive  flows.  However,  as  will  be  shown  in  this  section,  the  basic  lessons  that 
we  have  learnt  remain  important  when  it  conies  to  performing  simulations  of  detonation 
phenomena.  Thus  far,  for  simplicity,  we  have  utilized  the  so-called  Reactive  Euler  equations 
for  our  detonation  simulations;  a  single  reactant  A  is  converted  to  a  single  product  B  by  a 
one-step  irreversible  chemical  reaction  which  is  governed  by  Arrhenius  kinetics. 

In  one  space  dimension  the  reactive  Euler  equations  may  be  written  in  non-dimensional 
form  as 


(  P  ^ 

(  pu  \ 

(  0  \ 

d 

pu 

d 

p  u 2  -F  P 
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dt 
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+  3x 

( E  +  p)u 
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\  pz  1 

\  puZ  / 

V  KpZex p~h*fT  / 

Here  <>,  ti,  p,  T,  /?,  Z  and  E+  are  the  density,  velocity,  pressure,  temperature,  total  energy 
per  unit  volume,  reactant  mass  fraction,  and  the  activation  energy,  respectively.  Note  that  K 
is  a  free  parameter  that  simply  sets  the  spatial  and  temporal  scales.  Typically,  K  is  chosen 
such  that  for  a  ZND  wave  the  half-reaction  length  (the  distance  behind  the  detonation  front 
by  which  point  half  of  the  reactants  have  been  consumed)  is  scaled  to  unit  length.  The 
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following  equations  are  used  to  close  the  system  (2) 

E  =  pe  -f  pqZ  +  \pn2, 

V  =  [l~l)pe,  (2) 

T  =  p/p. 


Here  e  is  the  specific  internal  energy,  7  is  the  ratio  of  specific  heats  and  <7  is  the  heat  release 
parameter  for  the  chemical  reaction  A  — *  B. 

We  integrate  the  above  system  of  equations,  which  are  of  the  form  4^  -f  4-  =  s,  using 
the  method  of  fractional  steps 


w’,+i!  =  £sXc.£a.  w" 

The  source  operator,  £a,  corresponds  to  integrating  —  s,  which  in  this  case  reduces  to 

integrating  the  single  ODE,  42  -  -KZexpE+,T.  We  assume  that  the  temperature  field  is 
*  frozen  for  this  step,  allowing  us  to  use  the  nominally  exact  operator 

ZK+!  =  Zn  exp(  — A'  exp(  — E+/7’")Af), 

where  At  is  the  time  step  going  from  time  level  n  to  time  level  n  +  1.  For  the  convective 
operator,  Cc ,  which  corresponds  to  integrating  ~  -f  ||  =  0,  we  employ  Hancock’s  finite- 
volume  scheme  (Quirk,  1992a)  in  conjunction  with  the  adaptive  Riemann  solver  outlined 
in  this  paper.  Note  that  the  convective  operator  uses  a  time  step  that  is  twice  as  large  as 
that  used  by  the  source  operator.  The  generalization  of  this  integration  strategy  to  two 
space  dimensions  simply  consists  of  replacing  the  one-dimensional  convective  operator  by  its 
two-dimensional  counterpart. 

To  assess  the  capabilities  of  our  scheme  we  have  performed  simulations  of  one-dimensional 
pulsating  detonations,  or  so-called  galloping  instabilities,  for  which  bench-mark  results  ap¬ 
pear  in  the  literature  (Fickett  and  Wood,  1966;  Bourlioux  e t  al.,  1991).  Here  we  limit 
ourselves  to  presenting  the  case  where  the  overdrive  is  1.6  and  the  dimensionless  parameters 
that  appear  in  equations  (2)  and  (3)  are  given  by:  7  =  1.2,  £+  =  50,  q  =  50  and  K  =  230.75. 
We  have  run  this  problem  using  several  different  mesh  resolutions  in  combination  with  sev¬ 
eral  different  approaches  for  performing  the  reconstruction  step  of  our  Clodunov-type  scheme. 
Figure  7  shows  the  shock  pressure  history  for  the  case  where  the  computational  grid  pro¬ 
vided  40  mesh  cells  per  half-reaction  length,  and  the  reconstruction  process  operated  on  a 
characteristic  decomposition  of  the  conserved  variables.  Qualitatively,  this  pressure  trace 
is  identical  in  form  to  that  presented  by  Bourlioux  et  al.  (1991)  which  was  found  using  a 
scheme  that  fitted  rather  than  captured  the  detonation  front.  Figure  8  shows  a  convergence 
study  for  the  variation  of  the  peak  shock  pressure  with  mesh  spacing  for  two  popular  limiter 
functions,  namely  Superbee  and  Minmod  (Roe,  1985),  in  both  cases  the  reconstruction  pro¬ 
cess  operated  on  the  conserved  variables.  Note  that  a  relative  mesh  spacing  of  1  corresponds 
to  having  10  cells  per  half-reaction  length,  and  so  0.125  corresponds  to  having  80  cells  per 
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half-reaction  length.  Pleasingly,  the  grid -converged  value  for  the  peak  pressure  is  indepen¬ 
dent  of  the  reconstruction  process.  Moreover,  it  is  in  close  agreement  with  the  value  given 
by  Fickett  and  Wood  (1966),  and  the  ,alue  found  by  extrapolating  tin*  results  of  Hotirlioux 
ft  a/.(1991 ). 


Figure  7:  Shock  pressure  history  trace  for  a  calculation  with  40  mesh  cells  per  half-reaction 
length. 


Relative  Mesh  Spacing 


Figure  8:  Variation  of  peak  pressure  with  mesh  spacing. 

Even  for  this  relatively  simple  problem,  some  numerical  artifacts  can  occur  if  one  is  not 
careful.  For  example,  at  one  stage  in  the  cycle  of  the  pulsating  detonation,  the  profile  for  the 
reactant  mass  fraction  is  considerably  steeper  than  the  profile  for  the  steady  ZND  detonation 
wave  upon  which  the  mesh  spacing  is  based.  Consequently,  as  shown  by  Figure  9,  a  grid  which 
gives  10  mesh  points  for  the  half-reaction  length  of  the  ZND  wave  may  at  times  give  only  half 
the  expected  number  of  cells  for  the  pulsating  front.  For  a  compressive  limiter  function  such 
as  Superbee,  if  there  are  too  few  cells  covering  a  steep  but  continuous  profile,  the  profile  will 


appear  as  a  smeared  discontinuity  that  needs  steepening.  Hence  an  under- resolved  profile 
is  often  artificially  steepened  by  a  compressive  limiter.  Here,  such  oversteepening  causes 
anomalies  to  appear  in  the  shock  pressure  history,  see  Figure  10.  Note  that  such  anomalies 
are  not  associated  with  the  lead  shock  front,  and  so  the  question  of  whether  the  front  is 
fitted  or  captured  is  immaterial. 
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(b) 


Figure  9:  The  reaction  profile  for  a  pulsating  detonation  wave  may  at  times  be  considerably 
steeper  than  that  for  the  initial  steady  ZND  wave,  (a)  ZND  wave,  (b)  Pulsating  wave. 


Figure  10:  Anomalies  often  arise  when  an  under-resolved  calculation  employs  a  compressive 
limiter  function. 
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Following  Bourlioux  (1991)  we  have  also  conducted  a  two-dimen  sional  test  of  our  deto¬ 
nation  code,  namely,  the  simulation  of  the  transverse  waves,  or  so-called  cellular  structure, 
produced  by  a  detonation  wave  travelling  down  a  narrow  channel.  For  this  calculation  the 
overdrive  is  1.2,  the  channel  width  is  10  half-reaction  lengths,  and  the  dimensionless  param¬ 
eters  are  7  =  1.2,  q  =  50,  E+  —  10  and  A  =  3.124.  The  calculation  was  started  with  the 
solution  for  the  ZND  wave  as  initial  data.  This  planar  detonation  front  was  perturbed  by 
simply  allowing  it  to  ingest  a  small  region  of  fluid  where  the  rate  constant  l\  was  artificially 
decreased  by  20%.  The  calculation  was  then  run  until  the  ensuing  transverse  wave  structure 
was  fully  developed,  at  which  point  a  series  of  Schlieren-type  images  were  taken  so  that 
we  could  compare  our  results  with  Bourlioux’s.  These  Schlieren-type  images  are  shown  in 
Figure  11.  Note  that  the  calculation  applied  periodic  boundary  conditions  along  the  top 
and  the  bottom  of  the  channel  and  that  here  we  have  plotted  four  periods.  Qualitatively,  at 
least3,  these  results  compare  well  with  Bourlioux’s  results,  and  all  the  salient  features  of  the 
flow  have  been  resolved  (256  mesh  cells  covered  the  width  of  the  channel).  In  particular,  the 
regularity  of  the  repeated  vortex  patterns  are  the  same  for  both  sets  of  calculations. 

At  this  juncture,  it  is  worth  showing  the  results  that  were  produced  when  we  employed 
an  exact  Riemann  solver  rather  than  our  adaptive  Riemann  solver,  see  Figure  12  (only  one 
period  is  shown).  As  might  be  expected  from  having  seen  Figure  5  (d),  spurious  transverse 
waves  are  produced  in  the  vicinity  of  the  lead  shock  front.  These  spurious  waves  prevent 
the  correct  transverse  wave  structure  from  developing,  thus  ruining  the  simulation.  It  is  our 
contention  that  most,  if  not  all,  of  the  Riemann  solvers  that  are  commonly  employed  would 
be  similarly  unable  to  capture  the  lead  detonation  front  for  this  test  problem.  Whilst  some 
may  take  this  as  reason  enough  why  one  should  always  fit  the  lead  shock  front,  we  have 
shown  that  if  some  care  is  taken,  when  it  comes  to  detonation  simulations,  shock-capturing 
remains  a  viable  alternative  to  shock-fitting. 

In  an  attempt  to  unravel  some  of  the  dynamics  of  the  detonation  wave  for  this  problem, 
we  have  rerun  this  test  using  a  grid  that  was  four  times  finer  than  before.  In  order  to  achieve 
such  high  grid  resolutions  we  employ  a  relatively  sophisticated  mesh  adaption  scheme  the 
details  of  which  are  too  involved  to  give  here  (Quirk,  1991).  Figure  13  presents  a  pair  of 
Schlieren-type  images  for  the  temperature  field  which  show  the  local  flow  structure  just 
before,  and  just  after,  two  triple  points  collide.  It  would  appear  that  the  collision  gives  rise 
to  the  familiar  “explosion  within  an  explosion”,  which  in  this  case,  because  of  the  local  shock 
structure,  is  highly  anisotropic.  This  explosion  causes  slugs  of  hot  fluid  to  be  shot  fore  and 
aft  giving  rise  to  vortical  structures  which  are  similar  to  those  associated  with  Rayleigh- 
Taylor  instabilities.  Note  how  the  impact  of  the  forward  facing  jet  on  the  lead  shock  front 
causes  the  front  to  bulge. 

Obviously,  given  the  complexity  of  the  flow  field  for  this  test  problem  there  is  little  chance 
of  validating  every  detail  of  the  simulation.  However,  we  feel  that  our  computational  method 
has  matured  to  the  point  where  it  mav  be  relied  upon  to  provide  simulations  of  sufficient 
fidelity  for  fathoming  the  details  of  complicated  flow  mechanisms  as  is  done  here. 


3Owing  to  insufficient  information,  we  are  unable  to  perform  a  quantitative  comparison. 


Density 


Figure  Ii:  A  sequence  of  four  Schlieren-tvpe  snapshots  which  show  the  transverse  wave 
structure  of  the  detonation  front. 


Figure  12:  Spurious  waves  are  produced,  if  the  lead  shock  front  is  captured  using  an  exact 
Riemann  solver. 
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Figure  13:  A  pair  of  Schlieren-type  images  for  the  temperature  field  wliirli  show  the  sequence 
of  events  when  two  triple  points  collide. 


6.  Conclusions 

In  this  paper,  we  have  shown  that  Godunov-type  schemes  do  not  always  live  up  to  their 
reputation  as  being  models  of  robustness.  Specifically,  the  upwind  dissipation  provided  by 
almost  all  Riemann  solvers  is  inadequate  for  capturing  detonation  fronts  in  complex,  multi¬ 
dimensional  flows;  a  numerical  instability  can  develop  along  the  length  of  the  detonation 
front  which  interferes  with  the  genuine  transverse  waves  associated  with  the  propagation  of 
the  front.  At  the  very  least,  we  suggest  that  some  artificial  dissipation  mechanism  be  used  to 
augment  the  inherent  upwind  dissipation  so  as  to  suppress  this  type  of  numerical  instability, 
albeit  at  a  cost  of  some  loss  in  resolution.  In  general,  however,  to  improve  the  all-round 
robustness  of  a  Godunov-type  scheme  without  incurring  any  appreciable  loss  in  resolution, 
we  advocate  the  use  of  an  adaptive  Riemann  solver.  The  shortcomings  of  any  one  preferred 
Riemann  solver  are  circumvented  by  combining  it  with  one  or  more  complementary  solvers, 
such  that  an  individual  Riemann  solver  is  only  used  in  the  sorts  of  situations  for  which  it 
is  known  to  give  reliable  results.  Admittedly  this  approach  is  not  as  straightforward  to  im¬ 
plement  as  would  be  the  addition  of  an  auxiliary  dissipation  mechanism  to  one’s  favourite 
Riemann  solver,  but  it  does  prove  to  be  an  effective  means  for  producing  high  fidelity  sim¬ 
ulations  of  detonation  phenomena.  As  such,  it  provides  an  alternative  means  of  simulating 
detonation  phenomena  for  workers  who  might  otherwise  feel  compelled  to  fit  the  detonation 
front  solely  in  order  to  avoid  the  numerical  difficulties  associated  with  strong  shocks. 
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